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Abstract 

It has been established that there is an inherent limit to the accuracy of the reaction-diffusion 
master equation. Specifically, there exists a fundamental lower bound on the mesh size, below which 
the accuracy deteriorates as the mesh is refined further. In this paper we extend the standard 
reaction-diffusion master equation to allow molecules occupying neighboring voxels to react, in 
contrast to the traditional approach in which molecules react only when occupying the same voxel. 
We derive reaction rates, in two dimensions as well as three dimensions, to obtain an optimal 
match to the more fine-grained Smoluchowski model, and show in two numerical examples that 
the extended algorithm is accurate for a wide range of mesh sizes, allowing us to simulate systems 
intractable with the standard reaction-diffusion master equation. In addition, we show that for 
mesh sizes above the fundamental lower limit of the standard algorithm, the generalized algorithm 
reduces to the standard algorithm. We derive a lower limit for the generalized algorithm, which, 
in both two dimensions and three dimensions, is on the order of the reaction radius of a reacting 
pair of molecules. 
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I. INTRODUCTION 


Stochastic modeling has become a ubiquitous tool in the study of biochemical reaction 
networks [l-5|, as the traditional approach of deterministic modeling has been shown to be 
unsuitable for some systems where species are present in low copy numbers, or systems with 
spatial inhomogeneities [3, 6]. Instead stochastic, spatially homogenous or inhomogeneous, 
models are employed. 

Stochastic modeling can be carried out on multiple different scales. For processes occur¬ 
ring on the time scales typical of living cells we consider three different modeling scales: the 
spatially homogeneous well-mixed scale, the mesoscopic spatially heterogeneous scale, and 
the microscopic particle-tracking scale. In this paper the focus is on spatially heterogeneous 
modeling. 

A prevalent model on the mesoscopic scale is the standard reaction-diffusion master 
equation (RDME), in which diffusion of individual molecules is modeled by discrete jumps 
between voxels, while reactions occur with a given intensity once molecules occupy the same 
voxel. The next subvolnme method (NSM) 7] is an efficient algorithm for generating single 
trajectories of the system. The NSM has been implemented in several software pac kages, 


n 


including URDME [8|, PyURDME (www.pyurdme.org), STEPS (9|, and MesoRD [lO|. It is 
also available as a part of larger simulation frameworks such as StochSS (www.stochss.org) 
and E-Cell 11]. 

On the microscopic scale we model the molecules as hard spheres moving by normal 
diffusion. We track the continuous position of individual molecules, and molecules react 
with a probability upon collision. This model is commonly referred to as the Smoluchowski 
model [12], with the addition of a Robin boundary condition at the reaction radius of a pair 
of molecules. Algorithms aimed at accurately and efficiently simulating the Smoluchowski 
model for general systems have been implemented in E-Cell llj, Smoldyn 13], and MCell 

Q. 

It has previously been shown that there is an inherent bound of several reaction radii on 
the spatial accuracy of the RDME compared to the Smoluchowski model Q, Q . It was 
shown in 16] that by choosing correct mesoscopic reaction rates, the RDME could be made 
accurate all the way down to this lower bound. However, for mesh resolutions below this 
lower bound, the accuracy deteriorates. 
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In this paper we generalize the standard RDME by letting molecules occupying neigh¬ 


boring voxels react. Henceforth we refer to this generalization as the generalized RDME. 
Similar generalizations have been considered previously in 


cretizes the Doi model 


19] to obtain a convergent RDME. In 


5 , 


, 18|. In [17J, Isaacson dis- 
reaction rates are derived 


for a spherical model and applied to the RDME on a Cartesian mesh. In this paper we take 
a fundamentally different approach. By deriving reaction rates to match certain statistics 
of the Smoluchowski model, we arrive at analytical expressions for the reaction rates, and 
show that this approach yields accurate results down to a fundamental lower limit on the 
mesh size. This mesh size will be on the order of the reaction radius of two molecules. 


Importantly, we derive reaction rates under specific assumptions on the dynamics of 
dissociating molecules, and we show with a simple example that not doing so may lead to 
reaction rates that are inaccurate for certain systems. We thus argue that it is crucial to 
take dissociations into account in the derivation of reaction rates for the generalized RDME. 

The outline of the paper is as follows. In Section HI we review the Smoluchowski model 
and the RDME, and how they are connected through the mesoscopic reaction rates. In 
Section hh we describe the generalized algorithm, and derive accurate mesoscopic reaction 
rates as well as the lower limit on the mesh size. Finally, in Section IIVI we study two 
numerical examples, demonstrating the accuracy of the generalized RDME, and how it can 
be used to simulate systems that are intractable with the standard RDME. 


II. BACKGROUND 


A. Microscopic level 


At this level of modeling we track the continuous position of individual molecules moving 
by normal diffusion. Each species Si has a diffusion constant Di and a reaction radius a *. 
Consider two molecules, one of species Si and one of species S' 2 , with positions x ln and 

^ a 

x 2n at time t n . The molecules can react according to Si + S 2 W S 3 , where S 3 is some 

kd 

product. The probability distribution function (PDF) p(x 1; x 2 , f|x ln , x 2n , t n ) represents the 
probability that the positions of the molecules are given by x x and x 2 at time f; p then solves 
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the Smoluchowski equation 


d t p = -DiA Xl p + D 2 A X2 p. (1) 

It can be shown that with the change of variables 

Y = \J Lb-DiXi + a J D,D 2 x 2 (2) 

y = x 2 - xi, (3) 

we obtain two independent equations, where the equation for Y describes free diffusion, 
while the equation for y becomes 


d t p( y, t) = DA y p, 


( 4 ) 


where D = D\ + D 2 . Let a = (j\ + cr 2 be the sum of the reaction radii. We introduce a 
reactive Robin boundary condition at cr 

= k a p y (\\y\\ = a,t), (5) 

an l|y|| =cr 

where 


{ 2t raD, (2D) 
4na 2 D, (3D), 


( 6 ) 


and where k a is the microscopic reaction rate. The initial condition is given by p y (y, t n ) = 
6(y—y n ), and, since we assume that there is no outer boundary, we enforce p y (||y|| —* oo, t) = 


0. This equation can be solved analytically 20j, but the solution is difficult and expensive to 
evaluate numerically. Applying an operator split method to TPl - flS]) can significantly simplify 
the process of sampling new positions from the PDF 21]. 


An S 3 molecule is assumed to dissociate according to an exponential distribution with 
mean k d . Following a dissociation, the two products Si and S 2 are placed in contact a 
distance of a apart. 


A system of more than two molecules is not amenable to the direct approach of solving 
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for the full PDF, due to the high dimensionality of the problem. A common approach is 
instead to approximate the full problem as a set of one- and two-body problems, by dividing 
the system into subsets of single and pairs of molecules according to the distances between 
them. We can obtain a good approximation of the full problem by updating each subset 
independently during short time steps At. This algorithm is commonly referred to as Green’s 
function reaction dynamics (GFRD) [22, 23]. All microscale computations in this paper are 


carried out using a variant of the GFRD algorithm, described in 


21 ]- 


B. Reaction-diffusion master equation 


At the mesoscopic scale the simulation domain is discretized by N non-overlapping voxels, 
and diffusion is modeled as discrete jumps between the nodes of the voxels. The mesh may 
be either a Cartesian mesh, or an unstructured, tetrahedral (3D) or triangular (2D), mesh. 
A Cartesian mesh is suitable if the domain is simple, for instance a square or a cube, while an 
unstructured mesh has advantages for complicated domains. The jump coefficients between 
voxels are given by h 2 /(2dD) in the case of a Cartesian mesh, where h is the width of a 
voxel, d the dimension, and D the diffusion rate of the molecule. For an unstructured mesh, 
the jump coefficients can be obtained from a finite element discretization of the diffusion 


equation 


24], Reactions occur with some intensity when molecules occupy the same voxel. 


Let p(x, f|x n ,f n ) be the probability that the system is found in state x at time t, given 
that it was in state x„ at time t n . For brevity of notation, let p(x, t ) = p(x, t|x n , t n ). Let x*. 
and x. ; - denote the t-tli row and the j-th column of the K x S state matrix x, respectively, 
where S is the number of species of the system. The RDME is given by 


^ N M N M 

dt ■ P( ' X ’ ^ = EE °ir(Xj. - AOp(x 1-, ■ ■ ■ , X;. - fl ir , • • • , Xjy., t) ~ EE M X z>(x,f) 

i= 1 r=l 2—1 r= 1 

S N N 

+EEE djik (x.j ^ /, pfc)p(x.l, ••• , X.j Vijki • • • i X -Si t) 

j= 1 i= 1 k=l 
S N N 

-EEE ^ijfc(x.j)p(x, t), 

j= 1 i= 1 fc=l 

( 7 ) 
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where the propensity functions of the M chemical reactions are denoted by a ir (x t ), ^i ir are 
the stoichiometry vectors associated with the reactions, dijk are the jump coefficients, and 
Vijk are stoichiometry vectors for diffusion events. 

The RDME is in general too high-dimensional to be solved by direct approaches. An 
alternative approach is to generate individual trajectories of the system with stochastic 
simulations. The NSM [7] is an efficient algorithm frequently used for that purpose. 


C. Reaction rates for the standard reaction-diffusion master equation 


Consider a system of two molecules, one of species A and one of species B, that react 

ka 

according to A + B C, where k a and k d are the microscopic reaction rates. Assume 

kd 

that the molecules diffuse in a square (2D) or cube (3D) with periodic boundary conditions. 
Without loss of generality, assume that the B molecule is fixed at the origin, and that the 
A molecule diffuses freely with a diffusion rate D. The A molecule is initialized according 
to a uniform distribution. 

Let T meso (k™ eso , h) be the mean association time of the two molecules on the mesoscopic 
scale, and let T m i cro (k a ) be the mean association time on the microscopic scale. Under the 
assumption that r me so(^a ieso , h) = r m i CIO (k a ) holds, it is shown in Dll, 16] that the mesoscopic 
association rate is given by 


C“° = p [d] (K,h) = 


K 

h d 


1 + j j G«\h,a) 


-1 


( 8 ) 


where d is the dimension, 


G (d \h,a ) = 


£i°g(*-H)-j (i + A> ( 2D > 


jfe-S ( 3 D), 


( 9 ) 


and 


C d 


0.1951, d = 2 
1.5164, d = 3. 


( 10 ) 


The microscopic parameters are a, the sum of the reaction radii of the molecules, D, the sum 
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of the diffusion constants, and k ai the microscopic reaction rate. To simplify the notation 
somewhat, we let T^ eso (k a , h) := T meso (p ( ' d ' > (k a , h ), h). For a reversible reaction we match the 
mean binding time for h > h^, where 

{ ■v/vrexp a « 5.1a, (2D) 

V 4 t (n) 

|7rC 3 a ~ 3.2a, (3D). 

Let 'Tmeso ld (^a ies0 5 an d bScro d (^a) denote the average rebinding times—that is, the aver¬ 
age time until two molecules react, given that they have just dissociated—on the mesoscopic 
and microscopic scale, respectively. Again, to simplify notation, we let T™^ d,p (k a , h) : = 
r meso ld (p < ' d \k a ,h),h). The rebinding times can be written in terms of the average binding 
times 


T^o d ’ p (k a , h) = < eso (Ay, h) - r£ eso (oo, h) 

T mtro d ( fc a) = T micro ( k a ) ~ r micro (oo), 

where, for simplicty of notation, r^ eso (k a —» oo, h) and T m \ cro (k a —* oo) are denoted by 
r meso(°°;^) an d ' r micro(oo), respectively. That (TT2]l and (fT3l) hold can be realized by consid¬ 
ering the following argument. Given a uniform initial distribution, r meso (oo) is the time until 
the molecules are in the same voxel for the first time. By subtracting that time from the total 
binding time, we obtain the rebinding time. A similar argument holds for the microscopic 
case. We immediately see that because T^ eso (k a ,h) = r micro (fc a ) holds, the rebinding times 
will match if and only if t^^oo, h ) = r m i cro (oo). This holds for h = h^, and consequently 


( 12 ) 

(13) 


rebind, p/ 1 i\ . rebind (u \ 

'meso /£ 7 ^ 'micro 

for h > h ^ 

(14) 

rebind,p/ n U\ — _rebind/r, \ 
'meso \ rXj a-) ,L ) 'micro 

for h = h ^ 

(15) 

rebind,p /7 t\ . rebind /1 \ 

'meso \ 'micro 

for h < h*^. 

(16) 


As a mesoscopic dissociation event is a combination of microscopic dissociation and the 
diffusion required to get well mixed in a voxel, we require that should hold. 

For h < h ^ we cannot match the mean binding time while satisfying r ] !^ 1 “ d > T^ d “ d , and 
the accuracy of the RDME consequently deteriorates with decreasing h. Thus, h ^ is the 
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finest spatial resolution attainable with the standard RDME. 


For a given h > h^, we can compute the error in rebinding time as 


C“ nd,p (*:«, h) - = |w(oo, ft) - r mlc „(oo)|, 


( 17 ) 


where the right-hand side thus is a measure of how well-resolved a system is. Details of the 


above theory can be found in 


16], 


III. THE GENERALIZED REACTION-DIFFUSION MASTER EQUATION 


In the standard RDME, molecules react only when they occupy the same voxel. In this 
section we extend this approach by allowing molecules occupying neighboring voxels to react. 
To connect the standard RDME to the microscopic Smoluchowski model, we determined the 
rate with which molecules react when occupying the same voxel. For the generalized RDME 
we need to obtain the rates for molecules occupying the same voxel, but also the rates for 


molecules occupying neighboring voxels. In 15j, [16] we derive rates for the standard RDME 


by matching the mean association times on the two scales. To uniquely determine both of 
the rates for the generalized RDME, we need an additional constraint. 

In Sec. IIII Al we outline the algorithm, and in Sec. IIII Bl we derive mesoscopic param¬ 
eters by trying to match certain statistics of the microscopic model to the corresponding 
statistics on the mesoscopic scale. In Sec. IIII Cl we determine the dissociation rate of a 
reversibly reacting pair of molecules, and in Sec. IIII Dl we collect the results and summarize 
the algorithm. 


A. Generalized reactions 

Consider a domain D discretized by a Cartesian mesh, and a single reversible reaction 

^ ka 

A + B C. We extend the generalized RDME to allow reactions between molecules 
kd 

occupying neighboring voxels. Thus, if a molecule of species A occupies the same voxel as 
a molecule of species B , they react with an intensity given by k o. If the molecules instead 
occupy neighboring voxels they react with an intensity of k \, where two voxels are neighbors 
if they share one side. 














We can choose ko and k\ freely, with the restriction that the total intensity should be 
constant. Call the total intensity fc“ ieso . Let d be the dimension. Then, since each voxel has 
2d neighbors, k 0 and k\ must satisfy 


k 0 + 2dki = k™ so . 


Thus we can write 


k 0 = (1 - 2 dr)k™ so 
ki = rk™ so , 


(18) 


(19) 

( 20 ) 


where 0 < r < 1/(2 d). 

Now assume that a molecule of species C dissociates. We must determine where to place 
the two products A and B. It may seem natural to place them in the same voxel with 
probability 1 — 2 dr, and in neighboring voxels with probability 2 dr. While this arguably 
would yield the most accurate results compared to microscopic simulations for a single 
reversible reaction, we show below that this approach is unsuitable in general. 

First consider the single irreversible dissociation given by 

P ^ S 1 + S 2 . (21) 

In this case the microscopic and mesoscopic rates will be the same; thus /c™ eso = kd eg , and 
the products are placed in the same voxel. Now consider that in addition to (|2TT) we have 
the following reactions: 


sAsi ( 22 ) 

S* + S 2 S S 3 . (23) 

kd 

Again, (1221) is an irreversible unimolecular reaction and thus the mesoscopic and microscopic 
rates are the same. Now, if k* is large, the system f[2T]) - (l23j) will be well approximated by 

A' j 

P^*Sl + S 2 +± s 3 . (24) 

kd 
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Had we derived rates for reaction (1231) assuming that dissociating molecules are placed 
in neighboring voxels with some probability, we can see that the sequence (l24ll will be 
incorrectly simulated, as Si and S '2 are placed in the same voxel with probability 1 when 
P dissociates. Specifically, the rebinding dynamics of and S '2 will be incorrect, as the 
rebinding time will depend on whether they were produced from a dissociating S 3 or a 
dissociating P. 

To summarize: 

• Reactive molecules occupying the same voxel react with intensity (1 — 2 dr)k™ eso . 

• Reactive molecules in neighboring voxels react with intensity rk™ eso . 

• When a molecule dissociates, the products are placed in the same voxel with proba¬ 
bility 1 . 

The parameters r and /c“ eso now have to be determined from the microscopic parameters 
k a , a and D. 


B. Reaction rates 


Consider the reversible reaction A + BSC. Assume that the initial state of the system 

kd 

is given by one molecule of species A and one molecule of species B in a square (2D) 
or a cubic (3D) domain D of width L with periodic boundary conditions. For simplicity, 
and without loss of generality, assume that the A molecule is fixed at the origin while the 
B molecule has a uniform initial distribution and a diffusion rate D = Da + Db- On the 
microscopic scale the B molecule moves by continuous Brownian motion. On the mesoscopic 
scale, D is subdivided into non-overlapping squares or cubes of width h. The B molecule 
thus jumps between voxels with an intensity of kj = 2dD/h 2 , with d the dimension. Let 
'T’meso.r(fc“ ieso , h) denote the average time until the molecules react on the mesoscopic scale 
with the generalized RDME. 

For the local RDME, it was shown in [15, 16] that by enforcing the constraint r meso = r micro 
we obtain mesoscopic reaction rates as given by ([S]). I 11 addition, it was shown that 
approaches from above as h —> h Therefore it seems reasonable to require that, 

with the generalized RDME, we obtain an approximation of that is equal to or better 
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than the approximation we obtain with the standard RDME. The first constraint is therefore 
that the mean binding time agrees between the mesoscopic and the microscopic scales 


/ r meso j \ — 

l meso,r\ rx ' a i ,L ) 'ir 


)(^Q 


(25) 


and the second constraint will be that, given ([25]) . fc“ ieso and r minimizes the difference 
between the rebinding times at the mesoscopic and microscopic scales; that is, we want to 
minimize 


I rebind/1 meso U\ _ rebind/ t. \l 
|'meso,r', a >'"/ 'micro \ hj o-J \ i 


(26) 


under the assumption that (1251) holds, where r mes o, r is the average binding time (dependent 
on r and /c“ eso ), and where is the average rebinding time, in the generalized RDME. 

Note that with (1251) satisfied we have r meS o,o(^a Tieso ! = Aneso(&a ieso , h) an d T mesofi(^ es ° ■> M — 

rebind / Lmeso z, \ 

'meso \Aa ’i 11 )’ 


1. Mean mesoscopic binding time 

Again, assume that we have species A and B , with one molecule of each, and that the 
A molecule is fixed. The B molecule is initialized according to a uniform distribution, and 
diffuses with diffusion rate D. 

We start by deriving the mesoscopic mean binding time. To this end, let M* denote the 
average number of diffusive jumps required for the B molecule to reach a voxel at distance 
i from the A molecule, where the distance between two voxels is defined to be the smallest 
number of discrete jumps required to move from one voxel to the other. Let the set of all 
voxels at a distance i from the A molecule be denoted by di, let tj denote the average time 
for a diffusive jump, and let r* denote the average time for the B and the A molecule to 
react, given that the B molecule is occupying a voxel at distance i from the A molecule. 
Thus, tj = h 2 /{2dD) 1 and 


wff”. h ) = M ltj + T i- 


(27) 


The first term, Mjtj, represents the average time required for the B molecule to reach d\. 
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The second term, n, represents the remaining time until the molecules react, given that the 
B molecule occupies a voxel in d\. To obtain r mes0 , r we now derive analytical expressions for 
M] and T\. 


Lemma 1. Let N be the total number of voxels in the mesh. Then 


, f -K~ l N\ogN + {C'2 — l)N + 0(1) (2D) 

M] = { (28) 

[(C 3 -1)N + 0(VN) (3D), 


where C 2 and C 3 are defined in m- Let Ml denote the average number of steps required 
to diffuse from d t to dj. We have 


Ml 


N-2 
2d-l 


(29) 


Proof. First note that 


Ml = M° - M°. 


In 


26] it is shown that 


(30) 


n tt" 1 ^ log N + C 2 N + 0(1) (2D) 

M»=\ (31) 

[C 3 N + 0(y/N) (3D). 

Let be the average number of steps required to return to d 0 , given that we start in d 0 . 
The first jump of a molecule starting in d 0 always transfers the molecule to d\ , so we find 
that 


M° = M° - 1. 


(32) 


We know that Mg = N, shown in 25]. By combining (l30lh (l3l| . and (1321) . we obtain 


To obtain (129|) we note that we can write Mg as 


Mg 0 = 1 + — + 


1 

2d 


2d-l 


2d 


(Ml + Ml) 


(33) 
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To see that the above equality holds, start by considering a molecule in d 0 . The first jump 
transfers the molecule to dp, the second jump transfers it back to the origin with a probability 
of 1/ (2d), or to d 2 with a probability of (2 d — 1) / (2d). The average number of steps required 
to reach d 0 from d- 2 , is given by the average number of steps to reach d\ plus the average 
number of steps to reach d 0 , given that the molecule starts in d\. Now, solving (1551) for Mf 
yields ([29]). □ 


To obtain r meso ,r for < oo, it remains to determine T\. To that end, assume that the B 
molecule occupies a voxel in d\, and that the intensity with which the molecules react in d\ is 
given by l/(rfc“ ieso ). Then, to maintain a total intensity of l//c“ eso , the molecules must react 
with an intensity of 1/ [(1 — 2 d:r)k™ eso ] in do. We require that r > 0, and that 0 < 1—2 dr < 1. 
To simplify the notation we let l/(r/c“ ieso ) be denoted by pi, and 1/ [(1 — 2dr)fc™ eso ] by p 0 . 


Lemma 2. Let t\ be the average time until the molecules react, given that the B molecule 
occupies a voxel in di. Then 

(N P 2d- l)p 0 + (N + 2d- 2)t 7 - 
r i = - om -rTTT- -P 1 

2d(p 0 + t j ) +p i 
Po + tj N 
~ po + 2 drtj fc“ eso ' 

Proof. Let t° e and t\ denote the average time until the next event fires, given that the B 
molecule occupies a voxel in d 0 or d\, respectively. Then t° e = 1 / (pf 1 + tj 1 ) and t\ = 
l/ip^ + tj 1 ). 

By assumption, the B molecule initially occupies a voxel in d\. The next event can either 
be: (1) a diffusive jump, with probability pi/(pi+tj), or (2) a reaction event with probability 
tj/ipi + tj). 

Now assume that the next event is a diffusion event. Then: (1.1) the molecule jumps to 
d 2 with probability (2d — l)/2 d, or (1.2) the molecule jumps to d 0 with probability 1/(2 d). 
Assume that the molecule jumps to do- Then the next event is: (1.2.1) a reaction with 
probability tj/(po + tj), or (1.2.2) diffusion to d± with probability po/(po + tj). Thus, if the 
molecule is in the state (1.2), the time until the molecules react is given by 


(34) 

(35) 


ti 


To = 


< + 


Po 


Po + tj Po + tj 


(t° e + Ti) — t° e + 


Po 


Po + t 


-T\. 


(36) 
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Now instead assume that the molecule is in the state (1.1). The molecules cannot react until 
the B molecule reaches d\, and thus the average time until the molecules react is given by 


r 2 = M\tj + T\ 


N - 2 
2d — 1 


tj + Ti, 


(37) 


where M\ is given by (l29jh To summarize: 

The B molecule initially occupies a voxel in di, and the A molecule is hxed in do. 

(1) The B molecule diffuses with probability pi/{pi + tj). 

(1.1) The B molecule jumps to d 2 with probability (2d — l)/(2d). The average remain¬ 
ing time until the A and B molecules react is given by r 2 . 

(1.2) The B molecule jumps to d 0 with probability l/(2d). 

(1.2.1) The A and B molecules react with probability tj/(p 0 + tj). 

(1.2.2) The B molecule diffuses to d\ with probability po/(Po + tj). The average 
remaining time until the A and B molecules react is given by T\. 

(2) The A and B molecules react with probability tj/(pi + tj). 

Putting it all together, we obtain 


Tl 


tj | pi ( , 1 2d — 1 \ 

-—t e H-— ( t e + H-——t 2 ) . 

Pi + tj pi + tj y 2d 2d J 


(38) 


By inserting fl36l) and (1371) into fl38l) and solving for T\ we obtain ((34j) , after some cumbersome 
but straightforward algebra. By assuming N 1, (l35j) follows. □ 

Theorem 1 . Let r mes0jI . be the average time until the molecules react, given that the B 
molecules have a uniform initial distribution. Then 


fkr°,h) 


[n~ 1 N\ogN + (C 2 - 1)N] tj + 
(C3 - 1 )Nt, + 


N 


N 

t.meso 


po+4 rtj k% 

(3D). 


(2D) 


(39) 


Proof. This follows immediately from fl27jl . (|28jl and 


□ 


It is of interest to know the smallest voxel size h for which we can match the mesoscopic 
mean binding time, r mesor , with the microscopic mean binding time, r micro . In Q, 0 this 
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problem was solved in the case of the standard RDME for a general reversible reaction 

ka ... . _ t 

A+B ^ C. Similar results in the case of the generalized RDME can be obtained for the 

^d 

case of an irreversible reaction with k a — y oo. As we will find in Theorem 0 the lower bound 
for k a — y oo is in fact a fundamental lower bound for the generalized RDME. 


Theorem 2. Let h* ka be the smallest voxel size for which we can choose reaction rates such 
that r mes0!r (/C eso , h) = r micm (k a ). Then 


\ ^Fexp ( 3 + 2 AC 2 - iA qt ~ i,0599^ ( 2 D) 

Ko*={ V J (40) 

I f (^3 — l)vrcr « 1.0815a (3D). 


Proof. We do not have analytical results for r m i cro on a square or a cube, but given that 
L a is satisfied, an excellent approximation is provided by 


where 


,(K) 


j 1 ^ (2D) 

l£ ( 3D )- 


A = 7T2 


i a 


a = 


2vr D 


F{ A) 


log(l/A) 
(1 - A 2 ) 2 


3-A 2 
4(1 - A 2 ) ’ 


(41) 


(42) 


and where kcK = 4 naDka / (AircrD + k a ) is the classical mesoscopic reaction rate, valid for 
arg e volumes, derived by Collins and Kimball in 27]. The expression in 2D was derived in 


28], following the approach devised in [291. 


Since molecules are allowed to react with molecules occupying neighboring voxels, we 
obtain 


Ti —>• 0 for k a —> oo, r > 0 


(43) 
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and thus 


7meso,r M]tj for k a ->■ oo, r > 0. (44) 

We know M] from (l28lh and we have tj = h 2 /(2dD) by dehnition. We now obtain (1401) by 
solving 


M g tj — u m i cro (oo) 


(45) 


for h. 

In 3D, g5D becomes 


(C 3 - 1)L 3 


6Dh AttctD ’ 


(46) 


since M g ~ (C 3 — 1)N for N 1, and kcK —> 47t aD as k a —> oo. Solving (1461) for h yields 


h = -(C 3 — l)rca ~ 1.0815cr. 

o 


(47) 


In 2D, (l45j) becomes 


h?_ 
ID 


+(c2 -% 


= £ ny,2 

k a 2nD ’ 


(48) 


and for k a —> oo, we have 


k a 2 ttD 


fw L 2 , log (^)-l J . 2 


2ttD 


(49) 


In 


we used that A ~ 0 for I > cr. Now (1481) reduces to 


71 log I X ) + 


C, - 1 


= 7T 1 log I 7r 


iL 


cr J 47T 


(50) 


We can rewrite the equation above to get 


7r 1 log ^7T 2 


1 cr 

h 


l-c 2 3 
2 4tt' 


(51) 
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Solving for h yields 


3 + 27T ((f^2 — 1) 
4 


h = \/nexp 


^ a « 1.0599er 


(52) 

□ 


2. Mean mesoscopic rebinding time 


To satisfy the second constraint ()26l) we need both the microscopic_and the mesoscopic 
mean rebinding times. The microscopic rebinding time is derived in [h]], as 


rebind 
' micro 



(53) 


The mesoscopic rebinding time is simply given by 


rebind _ 

' meso ' 0 5 


(54) 


as To by definition is the time until an A and a B molecule react, given that they start in 
the same voxel. We have already derived r 0 in terms of T\ in (l36lh 


Theorem 3. Let be the average rebinding time of an A and a B molecule. Then 


rebind 
' meso,r 


~t° e + 


Po N 
p 0 + 2drtj kf eso ' 


(55) 


Proof. This follows immediately from Lemma [2] and (1361) . 


□ 


3. Solving for r and &:™ eso 

We now want r and fc™ eso to satisfy the constraints (1251) and (1261) . ft will prove useful to 
divide the problem into two cases: 

Case 1 : h > h ^ 

Case 2: h> 
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(56) 

(57) 




It turns out that in case 1 we get r = 0, effectively reducing the generalized algorithm to 
the standard algorithm. 


Theorem 4. Assume that r and k\ 
Then, for h > h*^, we have 


,meso 

a 


have been chosen to satisfy the first constraint (l25i) . 


rebind 
' meso,r 



(58) 


Since T™^ d,p (k a , h) > T^“ d (fc a ), it immediately follows that for h > h* OQ , the generalized 
RDME and the standard RDME agree. 

Proof. We already know that 



(59) 

(60) 


where r t , as previously defined, is the average time until the molecules react, given that the 
B molecule is in d % . The superscript g indicates that it is the average time in the case of the 
generalized RDME, and omission of the superscript indicates that it is the average time in 
the case of the standard RDME. 

We have assumed that (1251) is satisfied and consequently 


(M° - Ml)tj + (r 0 - rf) = Ntj + (r 0 - if), (61) 


0 = r T 


where the second equality follows from (1281) and (T3TT) . We know that = t 0 , so we get 



(62) 


Thus 


T r 


meso,r 


.rebind \ rebind 
meso.r — 'meso 


(63) 

(64) 

(65) 
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We have already shown that 


Now (165|) becomes 


Tn 


= t° e + 


Po 


Po + tj 


-T? > 


P 0 


PO + tj 


T? 


Po + tj N 

po + 2drtj k™ eso ' 


Po 


Po + tj N 

Po + 2drtj k™ eso po + tj p 0 + 2drtj k t 

Po + tj Po \ N 


p ° + ^ JL < N tj 

meso — J 
3 


Po + 2 drtj po + 2 drtj J k™ eso J 
t j N 

_ l _< ]\[t- 

po + 2drtj k™ eso ~ J 

1 1 

< 1. 


po + 2drtj k™ eso 


By definition, p 0 = 1/(1 — 2 dr)k™ eso , so (1711) becomes 


i 


(1—2 dr)k% 


2drtj k™ eso 


< 1 


meso 

a 


(1 — 2 dr)” 1 + 2 drtjk 
1 < (1 — 2dr )~ 1 + 2 drtjk 


< 1 


meso 
3 ,v a 


( 66 ) 

(67) 


( 68 ) 

(69) 

(70) 

(71) 


(72) 

(73) 

(74) 


Since 1 — 2 dr < 1, we have (1 — 2 dr) 1 > 1, and 2 drtjk™ eso > 0 so (1741) is satisfied for all r 
and £/f eso . Thus (1631) holds for all r and /c™ eso . □ 


What remains is to determine r and fc™ eso for h < h* 


Theorem 5. Assume that h = h* and that T\ t° e . Then 


/peso l\ 
1 meso,r \ A/ a 5 ' L ) 


T n 


o ( kd 


rebind / rmeso / \ 
'meso) ”7 


.rebind ft, \ 
micro 


(75) 

(76) 


for k™ eso = k a /h d and 1 — 2dr = 0. 


19 















For h < h 



rebind 
' meso,r 


(Jc meso M < H- reh[nd (k \ 

K^a -> ,L ) 'micro V^a)- 


,meso 


(77) 


Note that we have already shown that we can satisfy (125|) at least down to h = h*^ g . 
The assumption T\ t® means, in words, that the average time until two molecules react, 
given that they are one voxel apart, is much longer than the average time until the first 
event, given that they occupy the same voxel. Unless the microscopic reaction rate is 
very high, this should be a reasonable assumption for most systems. The necessity of this 
assumption is realized by considering two molecules in the same voxel. Now, if the average 
microscopic rebinding time is smaller than the average time until the first diffusion event 
on the mesoscopic scale, we could not hope to find mesoscopic rates that will yield a match 
between the mesoscopic rebinding time and the microscopic rebinding time. 

Proof. We already know that 


T, 


meso,r 



(78) 


and from assuming h = , it follows that 



(79) 


and thus 



(oo) + Ti. 


(80) 


Lemma 2 yields, for 1 — 2 dr = 0 and k\ 


,meso _ 

a 


k a /h d , 


N Nh d L d 


.rebind 

micro 


(*j. 


(81) 



We now have 



( 82 ) 
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and thus (1751) holds. Since we have assumed T\ t° e and 1 — 2 dr = 0, we get 


rebind 
' meso,r 




T 0 « T\ 


.rebind/r \ 
micro \™a)i 


( 83 ) 


and we have shown that (176]) holds. 

It remains to show (1771) . To this end, we simply note that 


Trr 


)() bre,so.r T M ' hnicro(^Ci) T T.) 


(84) 


since Mftj > T micro (oo) for h < . Thus 


rebind 
' micro 


‘Tmicro 


0«) 


'T’micro(oo) Tl 


_ rebind 

'0 • meso.r * 


(85) 


which proves (1771) . 


□ 


We now show that for h*^ < h < h ^ we can match both the mean binding time and 

the mean rebinding time. 


Theorem 6. Assume that h*^ g < h < h ^ and that n S> tg. Then we have 


/peso n 

/ meso,rv ri 'a "> ,L ) 


,(K) 


( 86 ) 


rebind / n meso l \ 
'meso ,rWa > ' V 


.rebind/n \ 
micro vAa/J 


(87) 


for 

k meso = f + fca/fr* \ K 

\tjQ 2 + Qk a /h d ) h d 
DQ(Q- 1) 

7 2 dDQ 2 + k a /h d ~ 2 ' 

where 


Q 


Ntj 

5ii m i' i( X I T meso r (oO 



(C 3 -1) 



(3D). 


( 88 ) 

(89) 


(90) 
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Proof. We already know that 


Trr 


r(C eso , h) = Mltj + Ti 


T mtro d ( fc a) = T micIO (k a ) ~ T micro (oo) = 


Nh d 


rebind / imeso 1 \ _ _ 

bneso,r ("'a i A) 


P 0 


AT 


p 0 + 2 drtj k™ eso 


Consequently we satisfy 


if 


M s tj T Ti — Tmicro(^a)) 


(91) 

(92) 

(93) 


(94) 


which, by (I34|) and (1351) . approximately holds if 


Po + tj N 
p 0 + 2drtj k™ eso 


Anicro ( k a ) M g tj 


To satisfy ([87)1 . we must, by (192)1 and (193)1 . satisfy 


Po N 


_ rebind/7 \ 

' micro V^aJ • 


(95) 


p„ + 2drtj fc”° 

Subtracting both the right-hand and left-hand side of (196)1 from (195)1 . we obtain 


(96) 


ti 


N 


Po + 2drtj k™ eso 


= r micro (k a ) - Mltj - T^Z d (k a ). 


(97) 


By definition, po — 1/(1 — 2dr)k™ eso and tj = h 2 /(2dD), so (197T) yields, after some straight¬ 
forward algebra, 


^,rneso _ 


(1 — 2dr)Ntj 


bnicro( k a ) - M]tj - T™ff(k a ) J rh 2 ( 1 - 2 dr) 


D 


(98) 


Since r mes0ir (oo) - Mltj and ) - 7 rili( . t , > (/oJ - T micro (co), (DU) becomes 


meso _ 


D 


Nti 


rh 2 V^micro(oo) - T me so,r(°o) 1 - 2 dr 


D 
rh 2 


Q 


l 


1 — 2 dr 


(99) 

( 100 ) 
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With fc“ eso as in (jlOOj) . we want to find r such that f|96jl is satisfied. Since = L d /k a = 

Nh d /k a , we obtain 


p 0 N Nh d 
p 0 + 2drtj k™ eso ~ If" 
1 1 _h d 

^ hdrtjpo 1 k™ so ~~ K 

<=> (1 + 2:rfrt jPo -')C~ = 

Since tj = h 2 /2dD and p 0 = 1/(1 — 2 dr)k a , (I103p becomes 


rh 2 

If 


(l-2dr)(C es °) 2 


+ k 


meso 

a 


t ,'>■ 


( 101 ) 

( 102 ) 

(103) 


(104) 


We expand the first term of the left-hand side to get 


rh 2 

IT 


(1 - 2dr)(fc” es °) 2 


D 

rh 2 


(1 - 2dr)Q 2 -2Q + j-W 


Thus 


(105) 


r h 2 n 

—(1 - 2 dr)(kzr°) 2 + *£”” = ^2 [(! - 2 dr)Q 2 - Q] , 

and (11041) becomes 


D 

rh 2 


[(1 - 2 dr)Q 2 -Q\ = 


h d ’ 


yielding 


DQ(Q - 1) 
2dDQ 2 + fc 


(106) 


(107) 


(108) 


Inserting r above into (HOOD yields (l88]h 

It remains to show that /c™ eso > 0 and 0 < r < 1/(2 d) hold for fc™ eso and r given by (188]) 
and (l89]h We first show that Q > 1, from which r > 0 follows. Thus we should show that 


Q 


NU 


> 1 


T, 


micro fOO) 7"meso,r fOO) 


(109) 
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holds. We start by showing that (11091) holds in 3D. By (I28)) . 


Trr 


r (oo, h) ~ (C 3 — l)Ntj = (C 3 — 1) 


L 3 h 2 

F 2 Id’ 


( 110 ) 


for N 1. We have already shown that 


L 3 


Tr 


licro(oo) — ' r meso,r(oO ) ^oo,g) ~ (^3 1) * \ 3 ^ji 

\^oo ,g) 


(ill) 


so (11091) becomes 


Q 


L 3 fe 2 
/ 1 3 2dD 


(ft -- (ft - i)| 


I, 3 Jf_ 
2 dD 


(Cb ~ 1) 

\ oo,9 


> 1 


(C 3 - 1) 


h 


k 


00 ,g 


-1 < 1 . 


( 112 ) 

(113) 


Since = C 3 /[C 3 — 1), and, by assumption, h < h^, we obtain 


(ft - 1) ( 7^- - l) < (ft ~V) f CS 


oo,g 


c, - 1 


1 = 1 . 


(114) 


Tims Q > 1, and, as a consequence, r > 0. In 2D we have 


Tmeso,r(oO, h) = [n 1 N log N + (C 2 - l)iV] tj 
T 2 T 2 T 2 

= 7T- 1 -£ r l0g^ + (C 2 -l) ^ 


/r 2 


2dD 


2dD ° h 2 


2dD 


(115) 

(116) 

(117) 


In 2D, similarly as in the 3D case, we have r micro (oo) = r mesO , r (oo, Kx>,g)- Thus 


Q = t 


LiJii- 

h 2 2dD 




L 2 


77 X 2dD lo §(7*^F 

1 




+ (C*2 l)2dD 2dD 1W 6 h' 

1 


lo sti + (ft - 1)535] 


2jr- 1 (log s #--logi) 2 ^ llo S 


(118) 

(119) 
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Since, by assumption, 


and 


we obtain 


1 < 


h 


h * 

n oo,g 


< 


hi 


h 


oo,9 


hi 


h 


= exp 


00,9 


3 + 2ttC 2 3 + 2n(C 2 - 1) 


71 

= exp [ - , , 


( 120 ) 


( 121 ) 


Q = 


> 


27t — 1 log 27r _1 log (exp f) 


= 1. 


( 122 ) 


Thus Q > 1 holds in both 2D and 3D, and we have r > 0. Note that with Q > 1, k^ eso > 0 
follows immediately. It remains to show r < 1/(2 d). To this end, we simply note that 


r = 


DQ 2 - DQ 
2 dDQ 2 + 3 


h d ~ 2 


DQ 2 - DQ 

DQ 2 +^I 2 d’ 


(123) 


where 


DQ 2 - DQ 

D Q 2 + 2ct^ 


< 1 


holds, since Q > 1. 

Thus 0 < r < 1/(2 d) and /c“ eso > 0 for h^ g < h < h* c 


(124) 


□ 


C. Dissociation rates 


Consider the same setup as before, with one A molecule and one B molecule reacting 

ka 

reversibly according to A+B C. Above we have determined how to choose the mesoscopic 

k-d 

association rates, so what remains is to determine the dissociation rate. This can be done 
completely analogously to the case of the standard RDME. We thus follow the approach of 
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161 ]. and conclude that we must have 


W°) 


-1 


T-rebind i (L.meso \ — 1 
'meso,r ' ) 



—rebind _i_ jC.—1 ’ 
' micro ' '' V:/ 


( 125 ) 


to obtain a steady state on the mesoscopic scale that matches the steady state of the micro¬ 
scopic scale. Thus it follows immediately that for /i^ < h < h we should have 


k 


meso 

d 


k d , 


(126) 


because r^(k™ s °, h) = r^ d (k a ) holds. 


D. Summary of the algorithm 


Assume that we have a cubic (3D) or square (2D) domain of width L , discretized by a 

ka 

Cartesian mesh with voxels of width h. Consider a reversible reaction A + B C, where 

kd 

k a and k d are the microscopic reaction rates. Let D = Da + D B , where Da and D B are the 
dissociation rates of species A and B , respectively. Let cr = a a + cs be the reaction radius 
of an A and a B molecule. 

The critical mesh sizes are given by 


{ 5.1a (2D) 

3.2a (3D), 


(127) 


for the standard RDME, and the critical mesh sizes for the generalized RDME are given by 


{ 1.06a (2D) 

1.08a (3D), 


(128) 


We now wish to simulate this system on the mesoscopic scale with the generalized RDME. 
The results of this section can be summarized as follows: 


For h > h*^-. The generalized RDME reduces to the standard RDME. Thus r = 0 
and fc™ eso = (k a , h). Molecules react only when occupying the same voxel. The 
dissociation rate is given by /c™ eso = h d k d k™ eso /k a , as shown in [lfj. 
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• For h^g < h < h*^. We match both the mean binding time, and the mean rebinding 
time, of the A and B molecules by choosing r and fc™ eso as in (I88|) and ([89]) . Now 
molecules react with an intensity of r/c™ eso when occupying neighboring voxels, and 
with an intensity of (1 — 2dr)k™ eso when occupying the same voxel. The dissociation 
rate is simply given by /c™ eso = k a . 

• For h < h ^ we can no longer match the mean rebinding time, and the accuracy 
deteriorates with decreasing h. 


IV. NUMERICAL RESULTS 


A. Rebinding-time distributions 


Consider a system of two species A and B , with one molecule of each. The A molecule 
is fixed at the origin, while the B molecule diffuses freely in space. In 0 it was shown 
that the local RDME matched the microscopic rebinding-time distribution for a reversibly 
reacting pair down to t* ~ (h*^) 2 / (2D). For t < t*, the behavior is inevitably going to be 
different, as the accuracy of the RDME is inherently limited by the spatial resolution. 

With the generalized RDME, we can match both the average binding times as well as 
the average rebinding times for h*^ > h > h* , and thus we could hope that also the error 
in distribution will be small at timescales of (h* oog ) 2 /(2D) < t < (h*^) 2 / (2D). 

In Fig. [j]we compare the microscopic rebinding-time distribution to the rebinding-time 
distribution for the generalized RDME. As we can see, there is a good match down to a 
spatial resolution of approximately cr. For the hnest meshes, the behavior at really short 
time scales is incorrect due to dissociating particles starting in the same voxel, but not 
reacting until they are in neighboring voxels. This introduces an error on the order of the 
voxel size, which will be on the order of the size of the molecules. 


B. Convergence of the generalized RDME 

The dynamics of some systems is resolved only at a fine spatial resolution. In particular, 
it has been shown that fast rebinding events can affect e.g. the response time of a MAPK 
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FIG. 1. In (a) and (b) we have plotted the rebinding-time distributions in 3D. For the standard 
RDME we have a good match between the microscopic and mesoscopic simulations for h ~ h*^, 
while the average rebinding time is underestimated for finer meshes. For the generalized RDME 
we see that the microscopic and mesoscopic distributions agree well for h^g < h < h ^ down 
to spatial resolution almost on the order of the size of the molecules, or a temporal resolution of 
approximately ( h* oo g ) 2 /{2D ). In (c) and (d) we have plotted the rebinding-time distributions in 2D. 
The conclusions are the same as for the 3D case. The parameters in (a) and (b) are given by a = 
2 x 10” 9 m, D = 2 x 10 -12 m 2 s -1 , L = 5.145 x 10” 7 m, and k a = 10” 18 m 3 s _1 . The parameters 
in (c) and (d) are given by a = 2 x 10” 9 m, D = 2 x 10” 14 m 2 s _1 , L = 5.2 x 10”' m, and k a = 
10” 12 m 2 s” 1 . Note that ‘standard RDME’ has been shortened to sRDME, and ‘generalized RDME’ 
to gRDME in the legends, to increase legibility. 


pathway 



We consider the system 


Jsi Hs 11 + s 12 ^s 2 

yS 2 -A 5*21 + S 22 -4- S 3 , 


(129) 
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which has a behavior similar to the MAPK pathway of jfi]. Due to the possibility of fast 
rebinding events, the long-term dynamics of the system is affected by spatial correlations 
between newly produced molecules. 


We start with an initial population of 100 Si molecules, with none of the other species 
present. The system is simulated for 2s, during which we sample the state of the system at 
201 evenly distributed points between t — 0 and t — 2. We simulate the system with both 
the standard RDME, as well as with the generalized extension, for different voxel sizes. Let 
S = {Si, Sn, S 12 , S 2 , S 2 1 , S' 22 , S 3 } . The error is then computed as 


m 


1 

201 


201 

EE II s ]” 

i= 1 ses 


[ 5 ] 


micro I 
i | 1 


(130) 


where [S]™' lcro is the average population of S at time fj, obtained with the microscopic GFRD 
algorithm, and where [S]^ so denotes the average population of S at time t; obtained at the 
mesoscopic scale with voxel size h. 

After a dissociation of either an S} or S 2 molecule from (II 29j) . the products can rebind 
quickly to produce an S 2 or S 3 molecule, respectively. On the microscopic scale, the products 
are in contact after a dissociation event, and thus the spatial correlation will be significant. 
At the mesoscopic scale, the products are placed in the same voxel after a dissociation. If 
the voxel size is large compared to the size of the molecules, the spatial correlation will be 
less than on the microscale. Thus, to simulate (112911 accurately, we would expect a fine mesh 
resolution to be required. 

Let a* be the reaction radius of molecule S',, and the reaction radius of molecule S t j. 
The parameters of the model are given by 


J k d = 10 s 1 
1 k a = 10 -19 m 3 s _1 


<j 1 = 10 -9 m 

(Til = <7 12 = 0.8 x 10 -9 111 
< a 2 = 2 x 10” 9 111 
C 21 = 022 = 1.8 x 10- 9 111 
u 3 = 2.5 x 10 ” 9 111 . 


(131) 


For simplicity, we let all species have the same diffusion rate, D — 10 12 m 2 s 1 . The S} 
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molecules are initialized uniformly in a cube of volume 10 18 m 3 . 

There is a critical lower bound on the mesh size associated with each of the system’s 
bimolecular reaction events 

[ hi : = /^(un + 012 ) ~ 5.0815 • 10 ~ 9 

{ (132) 

^h 2 := h^ 0 [o’ 2 1 + C 22 ) ~ 1.1433 • 10 8 . 

We know that for h > h* 00 , we are unable to match both the mesoscopic mean association 
time and the mesoscopic mean rebinding time to the corresponding microscopic quantities. 
Thus, for h > max{/ri, h 2 }, we will overestimate the rebinding time for both reactions, and 
consequently underestimate the average S 3 concentration. 

For h 2 > h > hi the dynamics is less obvious; we are underestimating the average 
rebinding time for the first reaction, but overestimating the average rebinding time for the 
second. As we can see in Fig. [ 2 ] (b), the positive and negative errors partly cancel out in 
this regime. At first the error decreases with decreasing h, but as we approach hi, it starts 
to increase again. The behavior of the standard RDME is hard to predict, and a priori we 
cannot be sure that a particular choice of h is suitable. 

In contrast, we see that the generalized RDME has a more predictable behavior, con¬ 
verging with decreasing h, and yielding an almost perfect match for h < hi. The difference 
in behavior is due to the generalized RDME matching the average rebinding time also for 
h < h*^, all the way down to h 2 , g : = h^ g {& 2 1 + 022 ) ~ 3.8934 • ICR 9 . 

V. SUMMARY 

For the standard RDME there is a lower bound on the mesh size, h below which the 
accuracy deteriorates. For h > h ^ we match the mean binding time of two molecules with 
the mesoscopic reaction rate given by p lyd \k ai h ). For h = h ^ we match both the mean 
binding time and the mean rebinding time of the two molecules. 

Some systems display fine-grained dynamics, requiring a fine spatial resolution to be 
simulated at the mesoscopic scale. By generalizing the standard RDME to allow reactions 
between molecules in neighboring voxels, we obtain a lower bound on the mesh size given 
by h^ g, where is on the order of the reaction radius of a pair of molecules. We derived 
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FIG. 2. In (a) we plot the average number of S 3 molecules as a function of time. As we can 
see, for a larger value of the voxel size h, we underestimate the number of S 3 molecules. For 
a very fine mesh, the number of S 3 molecules is overestimated. Somewhere in between we may 
obtain a good approximation compared to the microscopic results, but then the concentration of 
other species in the system will be incorrect. For simulations with the generalized algorithm, the 
average number of S 3 molecules is underestimated for coarse meshes, but as we refine the mesh, 
the dynamics approach that of the microscopic simulations. In (b) we see that the total error, as 
defined in (11301) . decreases down to hi for the generalized RDME, while error for the local RDME 
first decreases slightly but then increases as we refine the mesh further. We obtain an almost 
perfect match between the microscopic and the generalized RDME as h approaches hi, and all the 
way down to h- 2 , g . As in the legend of Fig. 1, we have shortened ‘standard RDME’ to sRDME, 
and ‘generalized RDME’ to gRDME in the legends. 


analytical expressions for the reaction rates, and showed that we match both the mean 
binding time and the mean rebinding time for h*^ < h < h*^. For h > h*^, the generalized 
RDME and the standard RDME agree. 

We studied the accuracy of the generalized RDME in two numerical examples. In the Erst 
example we showed that we not only match the mean rebinding time for h*^ < h < h*^, 

but that we also obtain a good match between the rebinding-time distributions at the two 
scales. In the second example we considered a system that cannot be accurately simulated 
with the standard RDME, as the mesh resolution required is below the fundamental lower 
limit h*^. We showed that with the generalized RDME we are able to simulate the system 
to high accuracy, and we showed how we obtain convergence to the microscopic simulations 
with decreasing mesh size h. 
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